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An approximate diagonalization method is proposed that combines exact diagonalization and 
perturbation expansion to calculate low energy eigenvalues and eigenfunctions of a Hamiltonian. 
The method involves deriving an effective Hamiltonian for each eigenvalue to be calculated, using 
perturbation expansion, and extracting the eigenvalue from the diagonalization of the effective 
Hamiltonian. The size of the effective Hamiltonian can be significantly smaller than that of the 
original Hamiltonian, hence the diagonalization can be done much faster. We compare the results of 
our method with those obtained using exact diagonalization and quantum Monte Carlo calculation 
for random problem instances with up to 128 qubits. 



I. INTRODUCTION 

Diagonalization of large Hermitian matrices is a diffi- 
cult problem in linear algebra, with applications in a va- 
riety of disciplines. In quantum mechanics, for example, 
the energy levels of a quantum system are obtained by 
diagonalization of the system's Hamiltonian. Knowing 
those energy levels is necessary for describing the behav- 
ior of the quantum system, e.g., the evolution of a quan- 
tum computer consisting of TV quantum bits (qubits). 
Calculating the exact energy spectrum of a Hamiltonian 
with usual numerical computation methods is possible 
only for up to TV ~ 20 qubits. For larger systems, the 
size of the Hilbert space (2^) becomes too large for the 
current level of available computer memory and speed. 

Although it is extremely difficult to calculate the exact 
spectrum of a multi-qubit system for large A'^, there are 
ways to calculate an approximate spectrum that shows 
important features of the exact spectrum reliably over 
a range of parameters. These methods include Density 
Matrix Renormalization [iHl] , Lanczos [5] , and Quantum 
Monte Carlo calculations [HHl]. Perturbation theory is 
also an approximate method, that is applicable when the 
system's Hamiltonian is close to a simpler Hamiltonian, 
called the unperturbed Hamiltonian, for which the eigen- 
values and eigenfunctions are known or easy to calculate. 
For example, the unperturbed Hamiltonian can be diag- 
onal in some known basis and the perturbation Hamilto- 
nian can have small off diagonal elements in such a basis. 
One can perform perturbation expansion in powers of a 
small parameter, characterizing the off-diagonal terms of 
the Hamiltonian, to find approximate solutions for the 
eigenvalues and eigenfunctions of the total Hamiltonian. 

A true perturbation expansion provides a Taylor ex- 
pansion of the energy levels in powers of a small param- 
eter. However, such expansion can become extremely 
complicated when there are energy level degeneracies in 
the spectrum of the unperturbed Hamiltonian. More- 
over, the perturbation expansion can quickly break down 
if the energy separations of the unperturbed states are 
small or when there are anticrossings between the eigen- 
states in the spectrum. Here, we combine perturba- 
tion expansion with exact diagonalization techniques to 



achieve an effective method for approximate diagonaliza- 
tion. The idea is to separate a subspace, e.g., low en- 
ergy states of the unperturbed Hamiltonian, from other 
(high energy) states in the Hilbert space. If the unper- 
turbed Hamiltonian is diagonal, then it might be easy 
to find its lowest energy states using the structure of 
the problem. Starting from the original Hamiltonian, 
we derive an effective Hamiltonian in the subspace using 
perturbation expansion. Perturbation theory brings into 
consideration, in the expansion of each term of the ef- 
fective Hamiltonian, the relevant states that are outside 
the subspace. If the unperturbed states in the subspace 
are non-degenerate, then there won't be a unique effec- 
tive Hamiltonian that can provide perturbed eigenvalues 
after the diagonalization. Instead, there will be an effec- 
tive Hamiltonian for each non-degenerate unperturbed 
state. Since the calculation involves exact diagonaliza- 
tion of these effective Hamiltonians, the final results are 
not true Taylor expansions in powers of the small param- 
eter, as in usual perturbation expansion. Yet, perturba- 
tion plays an important role in the derivation of these 
effective Hamiltonians. 



II. THE FORMALISM 

To derive the effective Hamiltonians, we adopt the pro- 
jection operator approach to perturbation theory as dis- 
cussed by Yao and Shi [10]. Consider the Hamiltonian 
H — Hq -\- V, in which Hq is the unperturbed Hamilto- 
nian and V is the perturbation Hamiltonian. The aim of 
the perturbation theory is to find the eigenvalues E„ and 
eigenvectors \n) of H, with 

iHo + V)\n) = E„\n), (1) 

assuming that the eigenvalues E^i^ and eigenvectors 
of Hq are known. 
Consider subspace S in the Hilbert space of Hq, con- 
sisting of Ng vectors \k'^^^). The subspace S can include 
both degenerate and non-degenerate eigenstates of Hq. 
We introduce projection operators 

P = ^ |fc("))(fc("'|, P = l-P. (2) 
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Since Hq, P, and P are diagonal in \k^^^) basis, we have 

[P,Ho] = [P,Ho]=0, (3) 

Let \k) denote an eigenstate of H that we are trying to 
calculate. We write \k) = \k)p + \k) p, where \k)p = P\k) 
and \k)p = P\k). We multiply both sides of (fTj), written 
for the eigenstate |fc), by P and P, respectively, to get 

El°^\k)p + PV {\k)p + \k) p) = Ek\k)p, 
Ho\k)p + PV{\k)p + \k)p)=Ek\k)p, (4) 

which can be rewritten as 

{Ek - 4"^ - PVP)\k)p = PVP\k)p, 
{Ek-Ho-PVP)\k)p^PVP\k)p. (5) 

Here, we have used P^ = P and P^ ~ P. Solving the 
second equation for |fc)p, we get 

\k)p^{Ek-Ho-PVPr^PVP\k)p (6) 

Substituting ^ into the first equation in ([s]), wc find 

H{k)\k)p^Ek\k)p, (7) 

where H is an x matrix defined by 

H{k) EE EfH + PVP 

+ PVP{Ek- Ho- PVPy^PVP, (8) 

with X being the Ng x iV^ unity matrix. 

So far, Eq. ([T]) is exact. What it means is that E^ is 
an eigenvalue and \k)p is an eigenvector of H{k) in S. 
Therefore, formally by diagonalizing H{k), one can find 
Ek and \k)p, but no information about other eigenvalues 
of H is obtained. Notice that Ek appears in both ([s]) and 
([t]) , therefore it has to be calculated self-consistently. As 
we shall see, perturbation theory can help to calculate 
Ek, order by order. 

If H is independent of k, then all the perturbed eigen- 
values and eigenvectors in S can be found in a single di- 
agonalization. For a /c-dependent H, on the other hand, 
only one of the eigenvalues, i.e., Ek, has physical mean- 
ing, while all other eigenvalues do not correspond to the 
correct eigenenergies of the spectrum of H. In that case, 
to calculate each Ek, one has to calculate its correspond- 
ing H{k), diagonalize it, and select the right eigenvalue 
that corresponds to Ek- Note that the states |fc)p found 
this way will not be orthogonal to each other. This indeed 
should be the case because only the original eigenstates 
|fc) in the full Hilbert space are supposed to be orthog- 
onal to each other, hence the projected states |fc)p may 
not be orthogonal. 

As mentioned earlier, H(k) is a function of Ek, which 
itself has to be calculated via diagonalization of H{k). 



The calculation becomes tractable using perturbation ex- 
pansion. Consider the expansion 

(Ek -Ho- PVP)-^ = (4°' -Ho- PVP + 6Ek)-' 

OO 

= E [{E^k^-Ho)-\PVP-SEk)]' {Ef^-Ho)-\ 

3=0 

where, 5Ek — Ek 



. Substituting into iSh , we get 



H{k) = E^^h + PVP + PVP X 



(9) 



OO 



^°'^-Ho)-^{PVP-5Ek)\ {Ef^-Ho)-^PVP. 
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Writing 5Ek = E\ 
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where the superscripts 



denote the order of perturbation, one can calculate H 
order by order. Because of the projection operator P, the 
operator P{E^^^ — Ho)~^P is not singular and is given 

by 



PiE^k^'-Ho) 



'P = 



E 



p(0) 



p(0)' 



(10) 



We now derive analytical formulas for all the elements 
of H{k) up to the forth order perturbation. Defining 
K,/3=(a(o)|V|/3(°)), where and denote un- 

perturbed states in S, and assuming that V has only 
off-diagonal elements in the chosen basis so that e'^'^ — 
(/c(o)|y|/c(°)) = 0, we find 
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Notice that for each added order of perturbation, a factor 
of the form Vnm/E^^^ is added to the expansion terms. 
The small parameter of the expansion, therefore, should 
be fc-dependent and have the form: 



A 



k ^ min [Vnm/El 

n.m4S 



(0)i 



(13) 



In our numerical calculations, we found the best agree- 
ment with exact diagonalization when expanding the di- 
agonal elements to the forth order of perturbation, but 
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the ofF-diagonal elements to the second order. This 
can be understood if one considers only two levels, i.e., 
Ng = 2. By diagonalizing the 2x2 reduced Hamiltonian 
corresponding to the two levels, if the diagonal elements 
are not the same, a second order off-diagonal element 
contributes a forth order term to the final eigenvalues, as 
it gets squared. Therefore to be consistent in the order of 
perturbation, one should expand the off-diagonal terms 
to the second order. 

Notice that all the energy differences in the denomi- 
nators are of the form and therefore depend on the 
unperturbed energy E'^"'' of state |fc^°)). As a result the 
calculated H is k dependent, unless all the states in S 
are degenerate (even in that case, the forth order cor- 
rection will still have fc-dependence through the second 
term in the last equation of (11)). As we mentioned be- 



fore, one cannot obtain all eigenstates by a single diag- 
onalization of H{k). Instead, one has to calculate H{k) 
for each unperturbed eigenstate \k^'^'^). The important 
task then is to select, among all the eigenvalues of H{k), 
the right eigenvalue that corresponds to state i.e., 
the perturbation of j/c*^")). If the perturbed levels do not 
cross each other, then Ek will be the fc-th eigenvalue after 
the diagonalization. In cases when the perturbed states 
do cross each other, the situation becomes more compli- 
cated. One can use the overlap of the new eigenfunctions 
with the old ones to identify which two correspond to 
each other. 

The accuracy of the calculations depends on the small 
parameter of the perturbation expansion, A^. To have 
an estimate of using (13), let E'min represent the 



lowest energy level outside the subspace S, therefore 
{Ekmlrn^s ^ Ei-ain—E^\ This providcs an upper bound 



(0)n 



for the small parameter: A/c < max^Vnm) / (Emin— E ^ 
Perturbation expansion, thus becomes more accurate for 
the lowest energy states for which E'^"' is smallest. Also, 
the accuracy of the calculations increases by increasing 
£"111111, i-fi-, increasing Ng. In principle, there is no limit to 
the accuracy and therefore no fixed radius of convergence 
as in the usual perturbation theory. By taking — )■ N, 
one can achieve 100% accuracy and an unlimited radius 
of convergence. In practice, however, Ng is limited by the 
limitation of the computation time and available mem- 
ory. By keeping TV^ small, the diagonalization can be 
done very quickly, but at the price of less accurate re- 
sults. Quite naturally, for < TV, small features of the 
spectrum that depend on the contribution of the higher 
energy states, beyond S and the states included pertur- 
batively, cannot be reproduced. 



III. COMPARISON WITH EXACT 
DIAGONALIZATION AND QUANTUM 
MONTE-CARLO SIMULATION 

To test our approximate diagonalization method, we 
study several Hamiltonians with different sizes and com- 
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FIG. 1: Hamiltonian parameters A and f as a function of 
normalized time s. 



pare our results with those of the exact diagonalization 
and quantum Monte Carlo simulation. The Hamiltonian 
we consider is an Ising Hamiltonian in a transverse field 
of the form 



H{s) 



-Hp = 



i<j 



(14) 



(15) 



where, s e [0,1], hi and Jij are dimensionless param- 
eters that can be adjusted, and A(s) and £{s) are en- 
ergy scales plotted in Fig. [l] Hamiltonian ( [m] ) was stud- 
ied in Ref. JL4, to investigate the scaling performance of 
an adiabatic quantum computation (AQC) [TT] proces- 
sor based on realistic Hamiltonian parameters. In that 
case, s = t/tf represents normalized time, where i/ is 
the total evolution time. It is known [ITj that in AQC, 
the minimum energy gap between the lowest two energy 
levels during the evolution determines the time of the 
computation. Therefore it is important to diagonalize 
the Hamiltonian ( 14 1 to calculate the minimum energy 
gap. Details of the AQC processor considered in this 
study are described in other publications [T?, 'T5'. Here, 
we only focus on the diagonalization of the Hamiltonians. 

In Ref. im a number of random Ising instances were 
generated and the size of the minimum gap in the spec- 
trum of their Hamiltonians was calculated using QMC 
simulation. With QMC simulation, one can calculate 
the energy gap between the lowest two energy levels us- 
ing the method discussed in OHl]. The instances used in 
Ref. [131 were generated by choosing hi uniform randomly 
from the set {-1/3,1/3} and a structured set of nonzero 
Jij values to be either -1, or uniform randomly from {- 
1/3,1/3}. The connectivity of the graph considered was 
motivated by a realistic quantum annealing processor as 
described in [T31 [H] ■ Here, we use the same set of prob- 
lems and compare our results with the QMC results of 
Ref.Hil 

First, we need to find the unperturbed eigenstates and 
eigenvalues of Hp for the instances studied. S ince l~Lp is 
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FIG. 2: Spectra of 8-qubit (a) and 16-qubit (b) random Ising 
Hamiltonians with transverse fields, relative to the ground 
state energy Eq. Solid thin lines represent the approximate 
diagonalization results and dashed lines represent exact diag- 
onalization results. Circles show the results of the quantum 
Monte Carlo simulations. 



already diagonal, it is only needed to determine the states 
with the lowest energy to form S. For that to be feasi- 
ble up to 128 qubits, we use the structure of the graph of 
nonzero couplings between the qubits. We use a dynamic 
programming method that is a variation on the bucket 
elimination algorithm [15 . We select an order in which to 
"eliminate" the qubits, i.e. to solve for the optimal value 
of a qubit conditional on the values of all qubits coupled 
to it that have not yet been eliminated. After eliminat- 
ing all qubits, the lowest energy state is simply retrieved 
by tracing back from the optimal value of the qubit that 
was eliminated last. To find the Ng low energy states, we 
keep track of the energy increase from choosing the sub- 
optimal value of each qubit being eliminated, and while 
tracing back through the elimination, the lowest energy 
Ns partial states encountered so far are kept, instead 
of just the optimal partial state. With this approach, 
the ability to find these unperturbed states is primar- 
ily limited by the largest number of qubits that need to 
be simultaneously considered during the elimination (the 
treewidth of the graph), instead of the total number of 
qubits. For the instances considered, the treewidth was 
up to 16. 
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FIG. 3: Approximate energy gap (solid lines) of example ran- 
dom Ising instances compared with results from Monte Carlo 
simulations (circles). From top to bottom N=32,48,72,96, and 
128. 



We choose the unperturbed and perturbation Hamil- 
tonians in the following way: 



(16) 



The small parameter in the perturbation expansion, 
therefore, is proportional to A{s)/£{s). In our calcu- 
lations we kept Ng ~ 2000 — 6000 in the subspace S. We 
choose iV^ in such a way that all the degenerate states in 
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the topmost energy level are included in the subspace S. 
A simplifying observation, for the calculation of (11), is 



that V only contains terms with operators of the form erf 
which flips the state of qubit i in the erf basis. For each 
V, therefore, only one bit flip from the original state is 
allowed. Consequently, in the calculation of the matrix 
elements of effective Hamiltonians H{k), using ( 11 ), only 
states with at most two bit flips from states |a)and |/3) 
participate in the sums. This significantly restricts the 
states |m), |n), or \p) that are summed over in (11). The 
states outside S were found by flipping qubits away from 
states in 5, as required by the perturbation expansion. 

We have examined many instances and compared the 
calculated approximate eigenvalues with the QMC simu- 
lation results of Ref. [T3] and also, for small size instances, 
with the exact diagonalization results. Here, however, we 
only report a few sample instances. Figure ^a.) and (b) 
show the calculated spectra for of 8-qubit and 16-qubit 
sample problems, respectively. Due to the small num- 
ber of qubits, the exact diagonalization was possible for 
these instances. The figures show excellent agreement 
between the exact (dashed lines) and approximate (solid 
lines) diagonalization methods over a wide range of the 
normalized time s. Even complicated details of the exact 
spectra are nicely reproduced by the approximate diago- 
nalization. As expected, some of the higher energy curves 
deviate from the exact diagonalization values at small s, 
where the perturbation expansion starts to fail. By in- 
creasing Ns, one can increase the validity range of the 
calculation at the expense of a longer computation time. 
The QMC results (symbols) for the above two instances 
are also plotted in the same figures. As can be seen, QMC 
agrees very well with the two other methods for s ^ 0.7. 
For larger values of s, QMC simulation fails to give re- 
liable results due to the small tunneling amplitudes. As 
we shall see below, the same pattern continues for larger 
scale problems. 

For problems with N > 16, it is not feasible to perform 
exact diagonalization as the size of the Hamiltonian be- 
comes exponentially large. As a consequence, we only 
compare our results with those calculated using QMC 
simulations. Figure |3] shows the calculated gap between 
the lowest two energy levels, for N from 32 to 128. For 
most instances the approximate diagonalization results 



agree with QMC calculations for 0.4 ^ s ^ 0.7. As be- 
fore, QMC fails to give the correct spectral gap for large 
s. Also, for the Ng values chosen, the perturbation ex- 
pansion becomes less reliable for s ^ 0.4, although the 
accuracy can always be enhanced by increasing Ng. In- 
teresting examples are = 48 and 128 for which there 
are anticrossings in the spectrum. The position of the 
anticrossing and the shape of the energy levels close to 
it are more or less consistent between the two methods 
of calculation. The size of the minimum gap at the an- 
ticrossing point, however, can depend on states outside 
the subspace S that are not included in the perturbation 
calculation. Therefore, the minimum gap size cannot be 
reliably predicted unless a very large number of states 
are included in S, or alternatively, the perturbation is 
expanded to high orders. 



IV. CONCLUSIONS 

We have developed an approximate diagonalization 
method for calculating low energy eigenvalues and eigen- 
functions of a large scale Hamiltonian. The method is 
based on derivation of a series of effective Hamiltonians, 
in a subspace consisting of low energy states of an unper- 
turbed Hamiltonian, using perturbation expansion. For 
each eigenvalue to be calculated, an effective Hamilto- 
nian is calculated and diagonalized separately. We have 
applied our method to find the energy eigenvalues of ran- 
dom Ising Hamiltonians in a transverse field. Our results 
agree very well with the exact diagonalization for 8 and 
16 qubit Hamiltonians, and with quantum Monte Carlo 
simulations for up to 128 qubits. The approximate di- 
agonalization method, however, is extremely faster than 
both of the above methods. 
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